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ON THE COMPENSATION OF GEOID ANOMALIES 


DUE TO SUBDUCTING SLABS 
ABSTRACT 

Candidate models of the forces which oppose the sinking of slabs are all constrained to 
produce results consistent with the following observation: relative geoid highs, which one assumes 
are due to slabs, characteristically occur over subduction zones. Certain models which are 
otherwise plausible, such as those based on a Newtonian half-space mantle, yield geoid lows instead 
of highs. This study has extended a published model of viscous corner flow in subduction zones 
in order to demonstrate that it can — in certain cases — produce the requisite geoid highs. Spe- 
cifically the relative geoid highs are produced if mantle flow is distinctly non-Newtonian (stress 
exponent n > 2). Results in the form of deflection of vertical (or geoid slope) profiles are com- 
puted for typical values of the slab parameters; they are compared with a representative profile 
of geoid slopes derived from Seasat altimeter data in order to show qualitative similarities. It is 
concluded that the effect of non-Newtonian flow as opposed to Newtonian, is to spread out the 
induced surface deformation, thereby stretching out the regional compensation to wavelengths, 
(transverse to the trench) of several thousand kilometers. 
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ON THE COMPENSATION OF GEOID ANOMALIES 


I 


DUE TO SUBDUCTING SLABS 


INTRODUCTION 

It is generally accepted that the lithosphere warms rather sluggishly as it descends into the 
mantle. As tlie cool slab is hydrostatically compressed it becomes relatively dense. A number of 

workers (McKenzie. 1969; Griggs, 197i; Kaula, 1972) have suggested that gravity (geoidal) highs 

observed over subduction zones may be due to the dense sinking slab. Recent work confirms 

this relationship between geoid high and subducting slab. Global scale modelling by Chase (1979) 

and Crough and Jurdy (1980) does so. A regional study by McAdoo (1981a) which uses Geos-3 

and Seasat altimeter data also supports such a relationship. 

Griggs (1973) stated and McKenzie (1969) implied that positive gravity anomalies (or geoid 

anomalies) due to the slab must be regionally compensated. Griggs (1973) also pointed out that the 
negative anomaly associated with the actual trench serves to partially compensate the positive 
anomaly due to the slab. In fact this partial compensation due to the trench can be expected to 
be more pronounced in a geoid height representation than in a gravity anomaly characterization 
such as tha^ used by Griggs. Trenches such as the Aleutian produce a short wavelength (X<400 km, 
transverse to trench) negative gravity anomaly; however they produce a corresponding negative 
geoid anomaly which has signifleant short and long wavelength (long wavelengths comparable to 
slab anomaly) components. According to McAdoo [1981b] the negative anomaly due to the 
Aleutian trench is not sufficient to fully compensate the geoid anomaly due to the Aleutian slab. 

In fact this study suggests that complete regional compensation involves induced surface deforma- 
tion with wavelengths transverse to the trench of several thousand kilometers. 

It is perhaps surprising that the broad geoidal highs typically observed over subduction zones 
do exist at all. More specifically, certain models suggest that the stress field induced by the sinking 
slab will give rise to a surface deformation and a concomitant geoid low which completely com- 
pensates the geoid high due to the relatively dense slab. Morgan’s (1965) model of a dense cylinder 
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sinking in a uniformly viscous mantle (Newtonian) produces, as a net result, a broad geoid low. Nu- 
merical experiments by McKenzie et al. (1980) on convection yielded geoid lows over the descend- 
ing flow limbs. Cold and hot blobs were introduced to drive the time dependent How. Davies' ( 1981 ) 
model of regional compensation of subducting slabs includes supporting stresses transmitted both up 
the slab and through the surrounding mantle. Just as the preceding models do, this model produces 
a net geoid low (s<n dashed line in Figure 4, Davies, 1981) for the case of complete compensation, 
i.e.. the case in which the magnitude of the supporting forces eipials that of the sinking forces 

A corner flow model like the one used by 1 ovish et al. (1978), Stevenson and l urner ( 1977) 
or McKenzie (1969) represents the kinematics of subduction more exactly than does either the 
sinking cylinder model (Morgan, 1965) or the convection/blob model (McKenzie et al., 1980). 

Tovish et al. ( 1 978) demonstrated that the induced flow in the mantle gives rise to dynamic 
pressures on the slab which may substantially counterbalance the gravitational torque on the slab. 
The gravitational torque alone would cause the slab lO hang straight down. Upon assuming that 
the torque due to induced flow exactly counterbalances the gravitational torque one can compute 
a net geoid signal. This paper is primarily devoted to computing such a signal. The computations 
indicate that a geoid low should occur over the subduction zones if the flow is nearly Newtonian. 
This result is in agreement with those of the two models cited above. However the computations 
also suggest that if the flow is substantially non-Newtonian a relative geoid high occurs over the 
subducting slab. The computed geoid highs have approximately the same wavelength as those 
actually observed over subducting slabs. 

THE CORNER FLOW MODEL 

As the slab founders it induces stresses in the surrounding mantle. These induced stresses 
give rise to a deformation of the surface which will affect the gravity field. Tovish et al. (1978) 
have developed a comer flow model of this induced flow. Their model will be used in this study 
to generate the induced stress field. The model provides analytical solutions for the induced stress 
field in a Newtonian or non-Newtonian mantle. More specifically, the appropriate constitutive 
equation is given by 
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t 


e = Tjn r- 


( 1 ) 


where e is strain rate, is a proportionality constant. : is stress and n is constant. By dclinition, 
n = I if the material behaviour is Newtonian. The model geometry and kinematics are shown 
in Figure 1. The slab is flat, its dip angle remains fixed, it is hinged (has infinite curvature) at its 
upper end, and subducts at a speed V^.. The mantle is assumed uniform in rheology , n are 
constant) throughout. Body forces internal to llte mantle which icsnlt liom li iii|>cialiiic vaiuii.His 
aie iicgicctcil as aie inertial loiccs. 

Tovish et al. obtain a solution for the induced stress field, Oj^, as follows. 


^^0 = 


P 

l^re 


^re 


where dynamic pressures p are given by 

p(r, 0) = TTo (2T?j,r)'^/" sin ((2n-l)'^ (d + 0„)/n)) 

and shear stresses by 

(2T?_r)'*/" 


T,e (r. 0) = 7T^ 


(2n-l) 


-n'/z 


cos ((2n-l)‘''-(0 + 0„)/n)). 


( 2 ) 


(3) 


(4) 


Constants it^ and are obtained by applying the appropriate boundary conditions (see Figure 1) 
to the velocity solution. Values of and 6^ are listed in Table 1 for different values of stress 
exponent, n. Actually these constants plus two others, 0 q and 0j, are necessary to completely 
specify the velocity field as given in Tovish et al. (1978). For completeness all four constants 
are shown in Table I. The constants are computed for n = 1 , 2, 3, 5 and a slab dip angle of a = 60°. 
They are obtained from four coupled nonlinear algebraic equations (excluding n = 1), and are 
solved numerically. Note that Tovish et al. also computed constants for n = 1 , 3 some of which 
are in error (G. Schubert, personal communication) and therefore do not agree with values in 
Table 1. Note also that certain of their equations (specifically A-IC, A-12, A-20 and A-21) contain 
minor typographical errors. The imaginary constants in Table 1 for the n = 2 case are easily 
understood if one recalls that constant tj^ in (1) becomes negative when shear stress, t. 
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Tab!.- 1 

Values of Constants in Corner Flow Stress and Velocity Field Expressions (dip a = 60°) 

Sec text for explanations. 




































is negative. When is negative and n - 2, the term in (3) and (4) becomes imaginary. 

For the case in which n = 2 and boundary conditions are coupled, the arc regime has a sector of 
positive shear stress and one of negative shear stress; therefore a pjecewise solution of the constants 
is required. 


SLOPE OF THE GEOID: ANALYTICAL EXPRESSIONS 

The corner flow model above does not include any stress or deformation conditions along trie 
lower boundary of the unsubducted lithosphere. In fact it will be assumed (as in McKenzie, 1969) 
that the lower surface (and upper surface) of the lithosphere deform such that the impinging normal 
stress vanishes. In other words the pressure equations, (2) and (3), determine the deformation of 
the lithosphere. The strength of the lithosphere will be ignored as the horizontal extent of the 
induced pressure field is substantially greater than the characteristic wavelength (~300 km) of the 
lithosphere’s response to transverse loads. Therefore the induced pressure p'(x) along the lower 
surface of the lithosphere causes a deformation, h(x), given by 

h(x) = p'(x, z = 0)/pg (5) 


where p is an average density for the lithosphere /upper mantle, g is the acceleration of gravity, 
h is positive upwards, p is positive in compression and x is defined in Figure 2 The 
gravitational effect of this deformation can be obtained by introducing an equivalent surface mass 
layer of mass per unit area. 


Op(x) = 


p'(x, z = 0) 
g 


( 6 ) 


The solution for induced pressure, p, as shown in (3) above may be substituted in (6) by using the 
following change of variables 


p'(x, 0) 


p(r = X, 0 = a), X > 0 
p(r = -X, 0 = ff - a), X < 0 . 


(7) 


It is assumed implicitly in equation (7) that slab thickness is 'effectively zero (cf. Figure 2). Note 
that for all cases covered in Table 1 pressure p'(x, 0) is everywhere negative (tensile) as is Op(x, 0). 
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The gravitational effect of the induced surface deformation in (5) or (6) can be obtained by 
evaluating the following integral expression for the slope of the geoid (or deflection of the vertical) 


dN 

dx 


(x, z = -d, 



Op(0 • (x-t) 

-E d^ 

(x-f)- + d- 


( 8 ) 


where G is the universal gravitation>:l constant, j.- is the iccderation of gravity and is the separa- 
tion distance between sea level and uppermost lithosphere as shown in Figure 2 It should be 
recalled that the corner flow mode! is two-dimensional (x, z) and extends to infinity in directions 
normal to the x, z plane. Therefore the surface mass layer exists everywhere on the z = 0 
plane and integrals having the form of (8) in many cases can not be evaluated (diverge). By succes- 
sive substitution of (3) into (7), tlien (7) into (6). and (6) into (8) one obtains a convergent integral 
for slope of the geoid (when n > 1 ). Upon evaluating the integral in (8) using equation 3.252. 1 2 in 
Gradsteyn and Ryzhik (1980) one obtains the expression for slope as follows; 


(x, z = -d) = d'* (x^ + d^) 

dx 


X cosec (7r(l - 1/n)) • 


Cj sin 


n 


where tj 


and where 
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s*” ((2n-l)^ (o:+ 0 ^, 3 )/n), 

^2 (7r-a + 0Qo)/n). 


(9) 
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Note that n^y^, 0^,^ and 0^^^^ arc the values of constants n^, 6^ in the expression for pressure (3) 
which apply in the case of the arc corner and oceanic corner respectively. 

The corresponding integral necessary for obtaining the geoid itself can not be evaluated. 

(i.e., it blows up). Similarly the integral in (8) diverges if d = 0. 

The expression (9), for geoid slope (denection of the vertical) is evaluated and plotted in 
Figures 3 and 4 for dip angle a = 60° and stress exponents n = 2, 3, 5. For the Newtonian case, 
n - 1 , the integral in (8) diverges and slopes can not be evaluated. Figure 3 represents the case of 
decoupled boundary conditions and Figure 4 represents the case of coupled boundary conditior.s. 

All curves in Figures 3 and 4 can be interpreted au broad (wavelengths on the order of 10^ km) 
geoidal lows. When represented as slopes (deflections of vertical) a geoid low characteristically lias 
the approximate shape of an inverted “N ” which has been streched horizontally and tilted counter- 
clockwise. This inverted “N” is, for example, evident in the n=2 curve of Figure 3. The curves in both 
Figures 3 and 4 have a gap between x=-250km and x= )-250km. In fact, expression (9) yields a continu- 
ous curve (the reader may connect the halves in his mind’s eye); however, the curves oscillate rather 
dramatically in this region (as do observed slopes; cf. Figure 8). Furthermore the model is unrealistic 
in this region due to the singularity in pressure, p(r = 0). and the assumption of an ignorable mechanical 
lithosphere. Values of the material constant and subduction rate, V^, are not specified to obtain 
the curves in Figures 3 and 4. instead the ratio is determined by postulating that torques 

on the slab due to induced pressure and gravity are in moment equilibrium. A uniform anomalous 
mass per unit area, o^, is specified for the slab as shown in Figure 2 and substituted into the moment 
equilibrium equation, 

-— a^gx^ cosa + I rpg (r, 0) dr 

^ Jo 

where p^ (r, 0) and p^lr, 0) are pressures given by equation (3) on the arc side and oceanic side, 
respectively, of the slab. In formulating ( 1 0) it is assumed that the thickness of the slab is effec- 
tively zero. By substituting (3) into (10) and rearranging one obtains the result, 
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Upon assuming that = 550 km /sin 60“ and that = 3 X 10^ gm cm'“ which is reasonably con- 
sistent with thermal model predictions (cf. McAdoo, 1981a, b, Schubert et al., 1975), the following 
values arc obtained for the sample case, n = 3 with decoupled boundary conditions: 


l/n 

I -= 1 .62 X 10* dyne cm'^ (12) 

wliere d = 6 km : 

Pg s 0.6 X 10^^ poises (13) 

where ii^ is effective viscosity defined as^i^ = and where t and are taken to be 

5X10^ dyne cm*^ and 5 cm/yr respectively. This viscosity is compared with those of other studies 
(e.g. Sleep, 1978) below. The quantity (Vj./2Tjj^d)’^" is needed to evaluate expression (9); it may 
be interpreted as a characteristic induced pressure 

Results in Figures 3 and 4 (also 6 and 7 below) are based on the assumed slab density Oj = 

3.0 X 10^ gm cm'^ ; however, results for any value of may be readily obtained by a simple change 
in vertical scaling because dN/dx, given in (9) is proportional to o^. 

Of course the geoid slopes given by (9) and c'lsplayed in Figures 3 and 4 do not include the 
effect of the anomalously dense slab itself. One can obtain the geoid slopes (deflections of vertical) 
due to the slab alone by evaluating the following expression 
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( 14 ) 


. , .-r’’ 

where sgn(x) = < 

I - I . X <0 

c = cosa . s = sina , 

and R“ = X" + r“ -2xr cosa. 

In the derivation of ( 14) the slab thickness is assumed to be negligible. The expression (14) is 
evaluated between limits r^ + 6r and 6r which are (cf. Figure 2) the distance from slab’s lower 
terminus to the trench and the distance from the slab’s upper end to the trench. Expression (14) 
is easily derived from equaMon (15) in Chapman ( 1979) by using a co-ordinate transformation and 
a derivative operation. Figure 5 shows an evaluated result from (14) using the same parameter values 
that were employed in Figures 3 and 4, i.e. a = 60°, o^ = 3 X 10^ gm cm'^, r^^ = (550/sin 60°) km, 

5r = (50/sin 60°) km. The curve in Figure 5 has the shape of stretched, upright, clockwise-tilted “N.” 
and therefore is representative of a relative geoid high. 

The composite of both effects, that of the slab and that of induced surface deformation may 
be obtained by adding the results of Figure 5 to those of Figures 3 and 4 respectively. These 
summations are shown in Figures 6 and 7. As in Figures 3 and 4 the results between x = -250 km 

and X =+250km are omitted for reasons discussed above; once again, the viewer may wish to 
connect the matching curve halves in his mind’s eye. For the nearly Newtonian cases (n = 2) 
one .--cs cither the equivalent of geoid low (stretched, inverted ‘N’ shape) in the decoupled 

case (Figure 6) or a relatively flat signal in the coupled case (Figure 7); this indicates complete 


9 


compensation. However for the n = 3 and n = 5 cases, particularly those with coupled boundary 
conditions (Figure 7), the resul:- are equivalent to a geoid high (i.e., they possess the characteristic 
stretched upright “N” shape). 

OBSERVATIONS AND DISCUSSIONS 

Geoid slope profiles obtained from model calculations can be compared with profiles of 
observed slopes. In general these observed profiles will be qualitatively similar to those model 
results based on assumed stress exponents n> 3 (strongly non-Newtonian). Figure 8 shows an 
observed sea surface slope profile derived from a pass ot Seasat altimeter data over the eastern 
Aleutian arc. The corresponding sea surface height profile is also included *n Figure 8. Ocean 
dynamics have not been removed with the exception of results from a global ocean tide model. 
However, most of tne long-wavelength (>500 km) content of the sea surface slope signal in Figure 8 
is geoid slope. The slope profile was derived using a runnmg, least square estimate of slope over 
four-second (5 point/27 km ) intervals. Ignoring the short-wavelength excursions in the observed 
slope curve ol' (Figure 8) one can see a stretched 'N' shape having the same sense, and roughly the 
same magnitude and wavelength as does, for example, the n = 3 model result in Figure 7. Large 
differences between absolute or average values of the model results (Figure 7;e.g, n = 5) and 
observed values (Figure 8) arc to be expected. The rather arbitrary choice of reference ellipsoid, 
mean versus hydrostatic, may strongly affect the average observed slope. Furthermore, the average 
value of model results is distinctly negative (e.g, n = 5). which is largely due to a model artifacts, i.e., 
infinite extent in x plus rectilinearity. 

Recall that the mode! calculations assumed that the slab dipped at an angle of 60°, extended 
to 600 km depth and possessed a density of = 3 X lO^gm/cm^. A dip angle of 60° is not incon- 
sistent with seismicity observations (see Jacob, 1972; Davies and House, 1979) for the Aleutian 
region covered by Seasat rev 553 in Figure 8; nor is the density estimate o^, inconsistent with 
inferred values (see McAdoo, 1981b; Grow, 1973) for this region. On the other hand, the slab 
model’s maximum depth of 600 km substantially exceeds estimates (see Jacob, 1972; Davies and 
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House, 1979) of the maximum depth for the Benioff zone in this region. It has however been 
suggested by McAdoo(1981b) that the anomalously dense slab may extend to greater depths than 
the scismically inferred 250 km. The implicit viscosity ~0.6 X 10~“ P. see above) is consistent 
with values used by Sleep (1978). In any case, theoretical and observational results show qualitative 
similarities; nothing more than this should be expected. 

On progressively longer spatial scales (e.g., those of order one earth radius) this infinite plane 
model provides an increasingly poor representation of the actual geometry and process involved. 
Therefcic a model based on spherical geometry is likely necessary to improve on the qualitative 
fits presented above. 

Examination of model outputs in Figures 6 and 7 reveals that large geoid lows flank the subduc- 
tion zone on the behind-arc side: somewhat smaller lows Hank on the ocean side. It is tempting to 
speculate that two major lows in the global geoid may in part he due to subduction-induced surface 
deformation. The first low extends from central Asia through India and the eastern Indian ocean 
(see Figure 2 or 3, Lerch et al., 1979), and generally lies a few thousand kilometers behind major 
island arcs of the western Pacific. The second geoid low exte.ids from the vicinity of the Bahamas 
southeastward across the general area of eastern Brazil, and lies behind the Middle American and 
South American arc. 

CONCLUSION 

The exact mechanism by which anomalously dense subducting slab' are supported has not 
been determined. However, candidate mechanisms must somehow preserve the characteristic 
relative geoid high which presumably results from the -.lab. Certain models based on a Newtonian 
half-space produce an induced surface deformation and concomitant geoid signals which completely 
compensate or wash out the relative geoid high due to the sinking slab. Tovishetal. (1978) presented 
a corner flow model which indicates that induced pressures may play a significant role in support- 
ing the slab. This study extends their model and demonstrates that characteristic geoid highs are 
preserved if induced flow is the primary slab support mechanism and if the mantle is quite 
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non-Newtonia,. (stress exponent n > 2). Non-Newtonian flow effectively spreads out the induced 
surface deformation, thereby retaining relative geoid highs and forcing part of the regional com- 
pensation process out to long wavelengths. The capability of the induced flow mechanism to yield 
a net geoid high strongly suggests that it is a primary operative process at subduction zones. 
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Figure 1 . Geometry, plate kinematics and boundary conditions for the 
corner flow model due toTovish ct al. (1978). Two cases are shown; 
coupled and decoupled* boundary condition. 



Figure 2. Geometry of the comer flow density model. 


15 




Figure 3. Geoid slope (or deflection of the vertical) versus transverse distance 
due to model of induced flow (alone). Decoupled boundary conditions are ass- 
umed. Three flow law cases (n= 2,3,5) included. 
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Figure 4. Geoid slope versus distance due to induced flow alone. Coupled 
boundary conditions are assumed. Three flow law cases included. 
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Figure 6. Geoid slope versus distance due to composite effect of induced flow 
(Figure 3) and slab (Figure 5). Decoupled boundary conditions assumed. 


19 








